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Abstract 

Holographic Particle Image Velocimetry (HPIV) is a novel technique for measuring 
the complete fluid flow around a body. Advances in computing power make this 
technique practicable for the first time. 

Currently popular techniques for experimentally determining fluid flow around 
a test body rely on measuring the flow at a single point and moving the sami)le 
point during the experiment. This implicitly time averages the data. HPIV can 
measure the flow field in a volume nearly instantaneously by sequentially “holo- 
graphing'’ the seeded flow field. 

Until recently the only practical method of processing the HPIV data was to 
physically reconstruct the optical image from the hologram and scan through the 
resulting volume using microscope optics. This is a very slow process. 

Recent advances in computing power make it possible to directly analyze the 
digitized holograph (digitized using a high resolution scanner) This reduces the 
post-processing time from several hours to less than 1-hour with much less expen- 
sive ecjuipment. 

The primary contribution of this thesis is an advanced multi-platform graphical 
user interface hologram synthesizer for use in calibrating digital HPIV reconstruc- 
tion algorithms. 

Thesis Supervisor: .lerome H. Milgram 
Title: Professor of Ocean Engineering 
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Chapter 1 
Introduction 

1.1 Overview 

This thesis presents an advanced iniiltiplatforni graphical user interface hologram 
synthesizer to produce calibration holograms of seeded flow fields. The synthe- 
sizer. HoloSynth, runs under X11R6 with OpenGL, and Microsoft Windows NT. 
Complete source code is contained in Appendix A. A brief introduction to parti- 
cle holography and holographic particle velocimetry is presented as well as results 
from computational reconstruction of the hologram generated by the synthesizer. 

1.2 Motivation 

Until very recently, particle image velocimetry (PIV) was limited to measuring 
fluid flow velocities only in a single plane. A sheet of laser light was passed 
through the flow which contained neutrally buoyant particles. A double exposure 
of the particles onto either film or a Charge Coupled Device (CCD) camera allows 
the two velocity components in the illuminated plane to be determined. 

The PIV method has recently been extended to three dimensions by Katz et. al. 
[1] by making a hologram of the flow field. Three dimensional Holograi)hic Particle 
Image Velocimetry (HPIV) uses either in-line or off-axis holographic techniques 
(see Section 2.2). Once the hologram has been recorded, tlie image is reconstructed 
l)y reversing the recording process and using a tliree dimensional optical seamier to 
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digitize the image for analysis by computer. This process can take several hours. 

An alternative to scanning the reconstructed image is to digitize the hologram 
and allow a computer to reconstruct the image in software. The recent advances in 
computing power should make this process far faster than the mechanical/optical 
scanning })rocess now used. 

Digitizing the hologram can be done in two ways. The initial method avoids 
the three dimensional scanning of the reconstructed image by using a high perfor- 
mance flat bed scanner to digitize an enlargement of the transparency produced 
in the normal holographic process. Alternately, a high re,solution scanner can be 
used directly on the hologram transparency. The flat bed scanning process is well 
understood, since the transparency can be enlarged to overcome resolution prob- 
lems with the scanner. However, it still requires recording the hologram on film 
with tlie attendant film processing. 

A more advanced technique would be to capture the hologram using a Charge 
Coupled Device (CCD) Camera. The hologram would be immediately available 
ill digital form, and no film or processing would be involved. 

Regardless of the hologram digitization technique the computational machinery 
for reconstructing the image must be calibrated. It far easier to produce a synthetic 
hologram for calibration than to set up a tightly controlled experiment and actually 
record a hologram. 

1.3 Requirements 

The initial goals of computationally reconstructed HPIV are to be able to accu- 
rately determine the flow field in a cubic volume a few inches on a side. The fluid 
speeds will be typical of those used in hydrodynamics research. The quantifiable 
requirements are: 

• Must be able to resolve a 10pm particle. This corresponds to the most 
conveniently sized neutrally buoyant seed particles used for Laser Doppler 
Anemometry. 
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• Accurate coverage within a 125cm^ (5ciii cube) volume. This will rise as the 
technicjne is perfected, and computing power increases. 

• Unambiguously detect fluid speeds < lOni/s. This accounts for free stream 
plus allowance for turbulence and model scale factors. 

The first and second recjuirements dictate the size and resolution of the recording 
device. 

Since holography is essentially a linear process, multiple exposures can be used 
instead of a rapid succession of exposures. So rather than design a system to 
rapidly ca[)ture two separate holograms, the recording medium is exposed two 
or more times without being altered between exposures. Tins results in multiple 
real images for each particle in the volume, the analysis process must identify tlie 
multiple images with a particular particle. This alleviates the need for a high 
speed film traversal mechanism and removes CCD output data rate as a concern. 

The third requirement sets the pulse lengths and pulse intervals of tlie laser 
used to generate the hologram. If the pulses are too long tlie diffraction patterns 
will be smeared from the particle motion. The greater the interval between pulses, 
the harder it will be to uniquely identify particle images with a single particle. The 
pulse length and interval requirements determine the recording device sensitivity 
and laser output power. 
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Chapter 2 

Review of Holographic Technique 

2.1 Overview 

A liolograiii is a two dimensional recording that can produce a three dimensional 
image. Typically electromagnetic (light) waves are recorded, but recording of 
acoustic waves and matter waves have also been made[2]. Electromagnetic holog- 
raphy differs from normal photography in that it encodes phase information as 
well as intensity information on the recording medimn. 

Recording media respond only to the intensity of the incoming wavefront, not 
its phase. To encode the phase using normal media, holography depends on using a 
coherent reference beam to interfere with a scattered or diffracted wave, or “object 
wave”, from the objects mider observation. As a result of this interference, the 
intensity at the recording plane depends on both the phase and amplitude of tire 
object wave. 

The original image is reproduced simply by illuminating the hologram with 
c-oherent light using geometry similar to that used to record the hologram. If the 
wavelength used for reconstruction differs from the wavelength used for recording, 
the image will be enlarged or reduced. 

The reconstruction process, while simple, is not amenable to rapid quantitative 
data collection. Quantitative analysis generally involves using three dimensional 
optical scanning of the reconstructed image. The ability to computationally re- 
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construct the image would dramatically increase the utility of holographic uiea- 
siiremeiit techniques . 

2.2 In-line (Gabor) Holography 

In-line holography, or Gabor Holography after its inventor, describes a holographic 
apparatus in which the reference beam is usually oriented perpendicular to the 
recording media and illuminates the objects from behind. It is not a useful tech- 
nique for objects that are large with respect to the recording medium, but if one 
is interested primarily in the three dimensional position of small objects and not 
their three dimensional shape, in-line holography can be very effective. It is com- 
monly used to size particulates and aerosols suspended in fluids, and it has been 
used for traditionally reconstructed HPIV[3]. 

Figure 2-1 shows the basic geometry and coordinate system for recording a 
plane in-line hologram. The object is illuminated by the reference source. The 
diffracted wave, or object wave, interferes with the wave from the reference source, 
or reference wave. This interference pattern is recorded by the medium in the 
recording plane and constitutes the hologram itself. 

If the reference and object complex wave amplitudes measured at the recording 
l)lane are given by R {x, y) and O (x, y) respectively, the intensity at the recording 
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plane is given by, 



I{x,y)= \R{x,y) + 0{x,y)\^ (2.1) 

Since the recording plane is perpendicular to the plane reference wave, the 
complex amplitude of the reference wave is a constant, r. Expanding (2.1) yields 

I{x,y) ^ \r + 0{x,y)f 

= r- + \0{x,y)\~ + rO{x,y) + rO* {x,y) (2.2) 

The intensity is recorded in any acceptable manner, usually as a positive trans- 
parency on film. Most reasonable recording media respond nearly linearly to in- 
cident wave intensity. The transmittance, t{x, y), of a positive transparency can 
be adequately described by 



t{x,y) = Tb-pTI{x,y) 



= n- pT 



+ \0 {x,y)\~ + rO {x, y) + rO* {x,y) (2.3) 



where is a constant background transmittance, T is the exposure time, and p 
is a constant of tlie recording medium. 

Figure 2-2 shows an in-line hologram being physically reconstructed. The 
wavefront transmitted by the hologram is simply the incident wave multii)lied by 
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the transmittance of the hologram and can be written, 



•0 (.r, y) = rt {x, y) 



= r {Tb + 0Tr‘^) + pTr\0{x,y)\~ 
+0T r^O {x, y) + 0T r^O* ( x, y) 



(2.4) 



The first term in (2.4) represents the nnifornily attenuated reference wave. The 
second term is normally very small, since typically O (x, y) r. It is called the 
"halo" and under proper conditions can be neglected. The third term is identical 
to, within a mnltii)licative constant, with the original scattered wave. It forms an 
image Zn from the hologram on the side opposite the observer. This is the virtual 
image, which is what the observer perceives. The fourth term represents a wave- 
front similar to the originally scattered wavefront, but with opposite curvature. It 
c'on verges to form a real image Zo from the hologram, between the hologram and 
the observer. [4] 

The observer then sees the virtual image, with a superimposed, oiit-of- focus 
real image along with the diminished reference beam. The two images instead of 
one and are the biggest drawbacks to in-line holograpliy. 



2.3 Off-axis (Leith-Upatneiks) Holography 



A simple way to separate the out-of-focus real image from the virtual image is to 
[)i'ovide a separate path for the reference beam and the object wave. Figure 2-3 
shows the basic geometry for recording an off-axis hologram. 

Because the reference beam is no longer perpendicular to the recording plane, 
it cannot be simplified to a constant, rather, it can be expressed as, 



R{x, y) ^ r exp{j27T^rx) 



(2.5) 



wliere G = (sin^)/A, since only the phase of the reference beam varies across the 
recording plane. 
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Figure 2-3: Off-axis liologram geometry 



111 a fashion similar to the derivation of (2.2) tlie intensity at the recording 
lilaiie is found to be, 

I (x,y) = \ 0{x,y)f + 2r |0(x,?/)| cos [27r^^x -h 0(x,y)] (2.6) 

Figure 2-4 shows the geometry for an off-axis hologram reconstruction. As- 
suniing the recording medium is linear, the transmitted wave complex amplitude 
is found in a similar manner, 

■0{x,y) = Ti,r exp{j27T^rX) 

+PTr \0{x,y)\^ exp{j2-K^rx) 

+(5Tr~0 (x', y) 

-\-(5Tr'0 {x,y)exp{jA'KS,rx) (2-7) 

Once again, the first term in (2.7) is simply the attenuated reference beam. 
The second term is a “halo” term, similar to the neglected term in (2.7). The 
third term is identical to the original scattered wave and forms a virtual image 
in the same position as the original object. The fourth term is, again, the real 
image, this time shifted from the optical axis by approximately twice the angle 
th(' rcfference beam makes with the recording plane. Hence the out of focus real 
image does not interfere with the vitual image. 
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Figure 2-4: Off-axis hologram reconstruction geometry 



At first blush it would seem off-axis holography is tlie answer to the virtual 
image problem for HPIV. If purely physical reconstruction is to be used, then off- 
axis holography is eminently practicable. The fatal flaw comes when considering 
digitizing the hologram for use in computational reconstruction. The reference 
beam now introduces a rapidly fluctuating intensity term that quadruples the 
resolution reciuirement of the recording medium and digitization system. For this 
reason off-axis holography will not be discussed further in this thesis. 

2.4 Diffraction 

Holography is fundamentally about diffraction. Goodnian[5] provides a detailed 
derivation of the Hnygens-Fresnel principle for a wave from constant 



recording plane. 

This expression, while physically meaningful, is not easily evaluated. Two 
fundamental simplifications are common depending on the geometry of the optical 
I)roblem being analyzed. The approximations are the Fresnel and the Fraunliofer 
aj)proxiinations. valid in the near and far fields, respectively. 




(2.8) 



where. 7’oi = + (x/, — x)^ -I- — yY, and Xh and are coordinates in the 



B 



I 



2.4. DIFFRACTION 



18 



2.4.1 The Fresnel Approximation 

A simpler and more useful expression of the Huygens-Fresnel integral can be found 
by simplifying roi using the paraxial api)roximation. . 

If 6 < 1. then, 



s/m> = \ + h- + 

Z o 



(2.9) 



Using (2.9) with b = ^ 



ro\ = \Jz^ ^ {xi, - x)- + {yh - y)"^ 




~ 2 






2 



(2.10) 



Substitution of (2.10) into (2.8) yields little simplification unless carefully cho- 
sen t erms are dropped. For the term in the denominator of (2.8) it is considered 
acceptable to drop all terms but 2 . However, small changes in the roi occurring 
in the numerator can cause significant changes in the value of the exponential 
and cannot be omitted. These large clianges are generated in part by the mul- 
tiplication of 7'oi by k, usually a very large quantity, and the sensitivity of the 
exponential to phase changes in its argument. The resulting simplification yields 
the Fresnel-Kirchoff Diffraction Integral: 



'0(^/1 , Vh) 




U (x, y) exp 




[(x;, - x)^ + {yh - yf] i dxdi2.U) 




— OO 






Where (2.12) differs from (2.11) only through factoring the term exp {x'f^ + yl)^ 
outside the integral. This approximation holds when the observer is in the near 
field, or Fresnel zone, of the diffracting aperture. Typical HPIV geometry will 
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place the recording medium in tlie far field of the individual particles but in the 
near field of the ensemble. 



2.4.2 The Fraunhofer Approximation 



Additional simplifications to the Fresnel Diffraction Integral can be made if the 
system meets the Fraunhofer approximation. 



2 » 



k{x‘^ -h y-)m^ 
2 



(2.13) 



Simplifying (2.12) yields. 






njkz 













(2.14) 



which is valid in the far field. The integral is recognizable as a two dimensional 
Fourier transform. 



2.4.3 Convolution interpretation of Fresnel-Kirchoff Inte- 
gral 

The Fresnel-Kirchoff Integral, 

k ^ 

ijh) = ^ jj U{x, y) exp [(x-/, - xf + {y^ - //)*] | dx dy (2.15) 

— OC' 



can easily be viewed as a convolution with kernel, 



h,{x,y) 



(,jkz 

jXz 



exp 






( 2 . 10 ) 



yielding 



fi’zU'hAUh) 



OC 

JJ U{x, y)h^{xh 

— OC 



^■,yh - y)dxdy = u{xhAyh) ® f>z{x, y) (2.17) 
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Tliis expression determines the field diffracted from a planar object, or an ensemble 
of co-planar planar objects. If a large number of very small iion-coplanar ol)jects 
are in the field of view. (2.17) can be extended: 



where the sum is taken over n planes within the volume. Each object plane. 

.(/). contains at least one diffracting object and is located z„ from the record- 
ing ])laiie. The planes need not be evenly spaced. 

2.5 Laser Characteristics 

2.5.1 Fundamentals 

At the atomic level, light is always generated by an atom transitioning from one 
energy state, Eq, to a lower energy state, Ei. When this transition occurs the atom 
emits a j)hoton with frequency n = (Eq — Ei)/h., where h is Planck’s constant. 
This emission is a probabilistic occurrence and may occur spontaneously provided 
the atom has sufficient internal energy. This is spontaneous emission. 

A second form of emission, first suggested by Einstein in 1917, is called stim- 
ulated emission. In the case of stimulated emission, a photon collides with an 
excited atom. The incident photon must have exactly the same energy as the 
transition the atom will undergo. When this condition is met, the atom will make 
the transition, emitting a second photon. The emitted photon is in phase with 
the incident photon and travels in the same direction. 

Stimulated emission is fundamental to laser operation. A material is induced to 
have a majority of its atoms in an excited state (called a “population inversion”). 
Since the majority of the atoms in the material are excited to the same state, 
they have identical transition energies. Initial spontaneous emission by a small 
number of the atoms will stimulate emission in others. The stimulated pliotons 
will in turn stimulate further emission in a run-away chain reaction. This in itself 
will not cause light amplification, but if the material is mounted within an optical 




(2.18) 



n 
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cavity (two parallel reflectors), the stimulated emission will preferentially orient 
itself along the axis of the cavity. Radiation not jiarallel to the axis will be lost 
to the surroundings and not contribute to the electromagnetic wave building up 
in the cavity. Since the stimulated emission is always in phase, the radiation 
traveling along the cavity axis constructively reinforces and the amplitude rajudly 
increases forming a longitudinal vibration mode witliin the cavity. If an end of 
tlu' cavity is only partially reflective, then some of the coherent amplified signal 
will be released, this is the usable laser beam. 

2.5.2 Establishing the Population Inversion 

The process of establishing the population inversion is called ‘pumping’. Pumping 
can be accomplished using flash lamps, usual in the case of ruby lasers, or electrical 
discharge, as in the case of most gas lasers. Photochemical reactions can also 
stimulate the population inversion but that technique is very rare. [6] 

The cavity material characteristics fundamentally determine what type of pump- 
ing will be appropriate. Continuous beam lasers must obviously use a sustainable 
pumping technique, such as tungsten filament incandescent lamps, or the afore- 
mentioned electrical discharge. 

The length of the pulse in non-switched pulse lasers is directly related to the 
length of the pulse length of the pump. Flash lamp pumped ruby lasers typically 
provide pulses on the order of several milliseconds long. The length of the pulse 
also varies from pulse to pulse as the condition of the flash lamp changes. 

The minimum time between lasers pulses in non-switclied lasers is directly 
related to the requirements of the pump. Flash lamps typically must cool imme- 
diately after use. The amount of cooling is directly related to the power output 
of the lamp, and can be on tlie order of minutes. This limits the pulse rej)etition 
freciuency of the unswitched laser. 

Switching can reduce pulse lengths to femtoseconds (1x10“^^ seconds) and 
pulse intervals similarly. Switching is discussed in some detail in Section 2.5.5. 
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2.5.3 Energy Levels 

When coiniuoii laser materials are pumped, the resulting atomic ]:>opulation usu- 
ally has many more than one po.ssible transition state available. The likelihood 
that an atom undergoes a particular transition depends on the quantuin mechan- 
ics of the particular material. In any event, when the population inversion starts 
to break down there can be a large number of longitudinal oscillation mod(«, 
('ach producing a separate wavelength in the cavity. Typically, a small number of 
closely .sjjaced wavelengths will predominate, but the resultant beam is not strictly 
inonochroinatic without s])ecial modifications to the laser (see Section 2.5.7). [6] 

2.5.4 Coherence Length 

The preceding discussion leaves the impression that all light generated in the laser 
cavity is coherent. This is not the case. Not all stimulated photons are descended 
from the same spontaneously emitted photon, and the Heisenberg uncertainty 
])rinciple places a lower limit on the precision with which the wavelengths are 
re])ioduced. So even a laser oscillating in only a single longitudinal mode has 
some degree of incoherence in its output. 

One measure of the quality of a laser’s output is termed “coherence length”. 
This length is the maximum beam path length within whicli the beam can coher- 
ently interfere with itself. The practical ramification of this is that the holographic 
system must have all path lengths from laser output to recording planes less than 
the coherence length of the laser. [7] 

2.5.5 Switching 

It is clear that the amplification properties of the laser cavity are critically depen- 
dent on the alignment and quality of the end mirrors. Misalignment or surface 
imperfections can destroy the amplification effect and turn the laser into an ex- 
])cnsive flash bull). 

This sensitivity can be turned to the laser users advantage. If the mirror prop- 
(n ties can be accurately controlled tlien the amplification process can be controlled 
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independently of the pumping process. 

Borrowing a term from electrical engineering, the ability of the cavity to res- 
onate is called it’s ‘Q’. Systems that switch the amplihcation process on and off 
by controlling the resonance qualities of the cavity are called ‘Q-switched’ or ‘Q- 
si)oiled’ lasers. 

Q-switching lias a snbth' benefit. The population inversion does not occur 
instantaneously at the beginning of the immping sequence. If the cavity cannot 
oscillate due to low Q, the population inversion can build to a much higher level 
than recinired for lasing. If the Q is suddenly increased, allowing the cavity to 
oscillate, the jieak energy in the pulse can be much higher than if the cavity was 
allowed to oscillate continuously through the pumping seciuence. So the jiower 
available in a Q-switched pulse can be much higher than in an nnswitched pulse. [6] 

An additional benefit of Q-switching is the ability to produce multiple pulses 
out of a single pumping process. This allows accurate control of continuously 
pumped lasers, and multiple rapid pulses out of short pulse pumped lasers such 
as ruby lasers. 

2. 5. 5.1 Electrooptical Q-Switching 

Several physical phenomena allow the electrical modulation of an optical signal. 
The most commonly used in Q-switching lasers is embodied in the Pockels Cell. 

Pockels Cells are formed from a crystal, typically potassium dihydrogen phos- 
phate (KDP). KDP has two distinct optical axes. With an electrical voltage 
applied across the crystal the index of refraction on one axis is different than the 
other. An incident polarized wave is split into two portions, one portion traveling 
more slowly than the other. The two waves are recombined upon exiting the cell, 
resulting in a wave that is elliptically, circularly, or linearly polarized de])ending 
on the applied voltage. 

The two voltages of interest for Pockels Cell Q-switching are termed the quarter- 
wavf'length (A/4) voltage, and the half- wavelength (A/2) voltage. The quarter- 
wavelength voltage will result in a circularly polarized beam at the output of the 
cell. The half-wavelength voltage will result in a wave linearly polarized 90‘” tf> 
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Figure 2-5: Pockels Cell Q-switch at (a) quarter-wave and (b) half wave retardation 
voltage. 



the input wave. 

The quarter-wavelength retardation system shown in Figure 2-5(a) is a ‘pulse 
off Q-switch. During the initial phase of pumping, the A/4 voltage is applied to 
the cell. The output of the Pockels cell is circularly polarized, after reflection the 
circularly polarized beam reenters the Pockels cell and emerges linearly jrolarized 
90° to the input beam. The polarizer then prevents the beam from reentering 
the laser rod, spoiling the amplification. When the voltage is removed, the Q is 
restored, and oscillation can start. 

The half-wavelength retardation system shown in Figure 2-5(b) is a ‘pulse on' 
Q-switch. With no voltage applied the Pockels cell does not alter the polariza- 
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Figure 2-6: Schematic Cavity Dump 



tion of the transmitted beam, and the crossed polarizers prevent oscillation from 
occurring. When switched on the Pockels cell shifts the polarization 90° allowing 
the beam to pass uninterrupted through the two polarizers. 

The half-wave retardation voltage for KDP is 7500V. KDP is also water soluble. 
Normally the crystal is hermetically sealed to prevent exposure to atmospheric 
moisture, but hydrophobic materials and extremely high voltages merit caution 
when used in a hydrodynamics experiment. 

2. 5. 5. 2 Mechanical Q-Switches 

There are two common methods of mechanically Q-switching a laser. They both 
rely on mechanically rotating a reflective element within the laser cavity. The most 
obvious is to rotate the mirror such that it is only periodically parallel to the fixed 
reflector at the oi)])osite end of the cavity. This system requires extremely critical 
alignment of the shaft upon which the mirror rotates. The alternative to rotating 
a mirror is to use a roof prism, which alleviates the alignment requirement. 

Mechanical switching is not as flexible as electrooptical switching, and is prob- 
ably not practicably capable of implementing the three pulse technique described 
in Section 3.1. 
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Figure 2-7: Cylindrical Transverse Electromagnetic Modes 



2. 5. 5. 3 Cavity Dumping 

Extremely short pulses can be generated using a “cavity duini)”. Cavity dumps 
allows the use of two 100% reflectivity end mirrors. When the laser j)ower is 
maximized, the cavity dump redirects the beam outside of the cavity, rapidly 
dumping the entire optical energy of the cavity into the output beam. The pulse 
length is reduced to the round trip transit time within the cavity. Cavity dumps 
can be used to switch intermittently and continuously pumped lasers. 

A schematic of a cavity dumped ruby laser is shown in Figure 2-G. Initially, no 
voltage is apirlied to the Pockels cell. When the flash lamp is fired the horizontally 
l)olarized wave is passed through the polarizer and not allowed to regenerate. 
Vertically i)olarized light is kept within the cavity and builds up the vibration mode 
witliin the cavity. When the energy stored in the Ruby has peaked, the Pockels 
cell is switched to its half-wave retardation voltage. The vertically polarized light 
is rotated to horizontally polarized light by the Pockels cell and passes through the 
Ijolarizer and out of the cavity. The combination of the polarizer, off-axis mirror 
and the Pockels cell essentially forms a voltage controlled mirror, whose reflectivity 
can be change from 100% to 0% and back to 100% during the cycle. [C, 8] 
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2.5.6 Transverse Electromagnetic Modes 

Laser cavities do not only oscillate longitudinally. There are also transverse inodes 
of vibration which are very similar to the vibration inodes of a drnni head. Laser 
specifications typically include a description of which transverse inodes are present 
in the output. The Transverse Electromagnetic Modes (TEM) are described using 
the symbol TEM,„„g. The first two indices describe the number of transverse 
inodes, while the third (often omitted in specification) describes the number of 
longitudinal modes. Figure 2-7 shows cylindrical transverse modes, along with the 
indices describing them. Rectangular modes are also possible. It is not difficult 
to find TEMoo ruliy lasers suitable for HPIV use. 

2.5.7 Output Frequency Control 

It is possible to dramatically reduce the number of longitudinal modes present in 
the outjnit. The simplest method is to use a Fabry-Perot etalon. This is simply 
a very small optical cavity that resonates at precisely the frequency desired in 
the laser cavity. The etalon effectively spoils the cavity Q for other frequencies, 
preventing them from building up. 

The simplest form of etalon is simply a thin glass plate. Internal reflections 
within tlie plate selection interfere with the wave to produce a sharp resonance 
peak. The effective thickness of the plate can be changed by rotating the etalon 
about an axis ])erpendicular to tlie cavity axis. 

Other etalon designs are possible. Alternate designs include wedge shaped 
plates that can be tuned by translating the plate, as well as designs that form 
^ul optical cavity using an air gap between two plates. Tlie latter design can 
sometimes be adjusted by varying the gas pressure or temperature within the 
cavity. [G. 8] 

2.5.8 Laser Construction 

Commercial lasers for field use are often manufactured in such a way as to limit the 
ability of the user to alter the laser. This increases laser reliability. Commercial 
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lasers for laboratory use are a different matter. As discussed above there are 
many additional components that can be addr*d to the laser cavity to customize 
and control the laser output. Laboratory lasers typically mount all components, 
including the mirrors and laser rod, on a long rail with additional space for adding 
polarizers, etalons. Q-switches and cavity dumps. 
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Chapter 3 

HPIV System Design 

3.1 Timing 

The requirements to unainl)iguously determine speeds less than lOm/s inside a 

volume 5 cm on a side impose severe restrictions: 

1. Exposure times must be short enough to ‘stop’ the motion. In a photograph 
the motion allowed diiring the exposure is related to the feature size of the 
object. During the photographic exposure the smallest features of an object 
should not move more than approximately 1/lOth their size or the image will 
be unacceptably blurred. Likewise an acceptable amount of particle motion 
during holographic exposure has been found to be approximately 1/lOth of 
a particle diameter [9]. 

2. The interval between exposures must be short enough to maintain the iden- 
tity of tlie particles, so tracking and velocity determination is possible. This 
restricts the mean free path length of a j)article between exposures to less 
tban the distance between particles. There is an additional problem of de- 
termining which particle image is ‘oldest’, so that direction of travel can be 
determined. 

3. The interval between exposures must be short enough to ensure a majority 
of the particles captured in the first exposure are within the covered volume 
during the following exposures. 
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For a 10/JiM particle, t liese reciuirements become: 

1. A 1/is exi^osure results in 1/lOth diameter of travel during the exposure. 

2. A lOps exposure interval results in 10 particle diameters of travel between 
exposures at lOin/s. 

Tlie ex])osure interval also places an upper limit on the concentration of seed 
particles. Assuring the particle mean free path between exposures is less than th(> 
average distance between {nirticles requires that 



Pp 



< 



|7r(100pm)"* 



2 X 10 



- particles 



cm-’ 



= 4 X 

in-^ 



Where p/, is the seed particle density, and the mean free path of the particles is 
ten particle diameters, 100pm. 

The direction i)roblem is subtler. Keeping the average inter-particle distance 
greater than the mean free path length between exposures signihcantly simplifies 
the task of identifying pairs of particle images as being the same particle. Unfortu- 
nately, even when the two images of a given particle have been identified, there is a 
180° ambiguity in the velocity of the particle. A simple solution to this j)roblem is 
to use three jjulses rather than two, with differing time intervals between pairs of 
pulses. For example, the second Ips pulse is sent 5ps after the first, then the third 
is sent lOps after the second. The pulse intervals are sufficiently short to permit 
j)article identification, while the difference in intervals allows determination of the 
t ime order of the images. An additional advantage of the triple exposure tech- 
ni(iue could be estimation of fluid acceleration. This technique assumes the flow 
speed will not change appreciably over the distances traveled by particle during 
the exposure sequence. If the particle were to radically decelerate to less than half 
of the sj)eed it traveled during the first set of exposures, it would be possible to 
confuse the time order of the exposures. 
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3.2 Recording and Digitization 

3.2.1 CCD Requirements 

The ultimate goal for computationally reconstructed HPIV is a completely digital 
process for determining tlie flow field throughout the volume. The most likely can- 
didate for the 'digital recording plane’ is a Charge Coupled Device (CCD) Array. 
Unfortunately, the resolution of CCD arrays is signihcantly less than holographic 
film. 

3. 2. 1.1 Characteristics of CCD arrays 

A CCD Array is a semiconductor architecture that transfers charge through stor- 
age areas. The CCD array is essentially an array of charge “buckets” . Photons 
are captured hy the buckets and converted to charge, the buckets only perform 
this conversion during a controllable integration period. More photons, or higher 
energy photons, hitting the array during the integration period result in a higher 
element charge. After integration the CCD Array support electronics scan through 
each bucket and report the charge level in each bucket. Tins charge is proportional 
to the intensity of the image focused on the array. [10] 

For holographic applications the fundamental characteristic of CCD Arrays 
that differentiates them from film is the limited number of recording elements. 
High quality photographic film can have grains smaller than Ipni over an area 12- 
20 centimeters stjuare, or better than 120.000x120,000 recording elements. The 
best CCD arrays currently commercially available have element sizes on the order 
of 10 // 111 , but the total array size is only about 5 centimeters on a side, yielding 
an array of approximately 5000x5000 elements. This is tlie root of the limitation 
in holographic applications. 

3.2. 1.2 CCD Resolution Requirements 

3. 2. 1.2.1 Diffraction patterns of circular objects Parrent and Thomp- 
son analytically derived the diffraction pattern for an opaque circular object in a 
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coherent, field. [11] For a circular particle, the field and intensity, respectively, are: 



Where, r is measured from the center of the object projected into the recording 
l)lane. and A(q) = .7i(tt)/a, a Bessel function of the first kind divided by its 
argument. 

An illustrative example is plotted in Figures 3-1 and 3-2. The parameters 
more characteristic of HPIV systems result in plots that do not acceptably show 
the salient features of the intensity distribution. 

3. 2 . 1.2. 2 Recording the Diffraction Pattern Equation (3.1) describes the 
diffraction pattern recorded for a single particle with circular cross section. All 
previous discussion has assumed this diffraction pattern is perfectly recorded. In 
rc'alitv the recording medium has limited bandwidth (he. a grain), and limited 
physical extent. 

Figure 3-1 shows a high fre(juency sine function modulated by a lower frequency 
A(a) envelope function. The envelope function is dependent on both particle 
diameter and particle distance from the recording plane. The high freciuency 
sine function is dependent only on particle distance. Accurately recording the 
high frequency component is the most critical factor in obtaining good depth 
discrimination. 

Vikram[3] recommends recording the first three side lobes formed by the bessel 
envelope function. The spacing of the fine fringes. Ax, in the third lobe can be 
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Figure 3-1: Ditfractioii iiiteiisity of a rirnilar ohjert with a — 50//in, A = ().G913//m, and r ^ 10( in 
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Figure 3-2: Enlarged plot, of the third side lobe in Figure 3-1 



found using the arguinent of the sine term in (3.1), kr^ j2z\ 



kx\ kx\ 

2 2 
Xo — x\ 

Ax 



2tt 

2xAx — 2Xz 
Xz 



(3.2) 



Equation (3.2) should be evaluated in the region containing the third side lobe. 
The mid-radius of the third side lobe, is: 



•A i^n) / 

C\ji 

k'n 
rz 

Equations (3.3) and (3.2) give fs = 0.189cm and Ax = 3.65 x 10“^cm for 
the iiarameters plotted in Figures 3-1 and 3-2. Using 2 = 5cm and a - 10//m, 
with A = 690mn. more reasonable for HPIV purposes, 7~3 = 0.472cm and Ax — 
7.31 X KrU-in. 

The parameter dictates the minimum size of the hologram, with no mag- 
nification. and is easily met. The parameter Ax dictates the minimum grain of 
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the hologram recording medium. Since enlargement or reduction of the hologram 
during recording is easily accomplished using lenses, the absolute size of the fringe 
spacing is not important. Rather, the important quantity is the ratio the fringe 
si>acing to the total size of the system, so that the fine fringes produc:ed by i:>ar- 
ticles at opposite extremes of the volume can be resoh'-ed. The sampling theorem 
suggests that the recording medium should be able to sample at twice tlie high- 
est spatial frequency required, or 1/2 the fringe spacing. To achieve the desired 
accuracy, a CCD array would have to ha\'^(’ 

2 • 0.05cm • 1 

; ~ 13680 i)ixels square. 

7.31 X lO-h-m 

This resolution will adequately record the particle diffraction patterns. Unfortu- 
nately no CCD arrays are commercially av'-ailable that reach this resolution. The 
best currently available are in the range of 5000 pixels square. 

3.2.2 Exposure Requirements 

While the performance characteristics of CCD arrays are well known it is difficult 
t o j)i'edict the amount of power incident on tlie array from the suspended partic'les. 

The resolution analysis was performed assuming the finest necessary fringes 
were detected. The incident power depends heavily on the power output of the 
laser, the absorption of the fluid, and the reflectivity of the particles, and the 
viewing window material. 

The exposure problem is different in character from tlie resolution problem. 
The resolution problem has to wait for better CCD technology, or lower the ac- 
curacy reciuircments for HPIV significantly. The exposure problem can be solved 
many ways. the. simplest (but not necessarily the least expensive) way is to increase 
the laser power. 
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Chapter 4 
Simulation 

4.1 Motivation 

A primary factor preventing HPIV from becoming an easily and widely used tecli- 
iiicpie is tlie data reduction effort. As previously mentioned, physical reconstruc- 
tion techniques are simple, but require costly equipment, and inordinate amounts 
of time. If computational reconstruction methods can be perfected, HPIV could 
be Ijecome more practical. 

There are several possibilities for computational reconstruction algorithms, 
one of which (Onurars method) is presented in Section 5.1. The methods must 
be evaluated with realistic input in order to gauge their effectiveness. Obtaining 
realistic input is difficult. There have been experiments done using holographic 
techiiiciues to study small ocean biological specimens, but obtaining the actual 
holograms produced is difficult. An additional difficulty is that you don’t know 
precisely what was contained in the original image so calibrating the reconstruction 
algorithm is nearly impossible. 

4.2 Hologram Synthesis 

There are several methods applicable to synthesizing holograms similar to those 
exj^ected from HPIV. Tlu> most straightforward methods calculate the individual 
r'ontribution from each particle to each pixel on the simulated recording plane. 
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Tlie.se inetliods are the mo.st accurate and flexible but suffer from excessive com- 
putational reriuirements when large numbers of particles are lieing simulated. 
One way around this problem is hinted at liy Ecpiation (2.18): 

y) = '^Un (.X, y) (8> (x, y) 

n 

W'e can subdivide the region being simulated into a large number of thin slabs. 
Within each slab there can be a large number of particles, but each will be treated 
as if it had the same coordinate. We tlien calculate the contribution of each 
slab to the total hologram rather than the contribution from each particle. With 
a larg(' number of jiarticles we can achieve significant performance increa.ses with 
little degradation of quality. 

Evaluating (2.18) is quite efficient. Using the convolution theorem: 



■0 (x, y) 



^ Un (x, y) (g) (x, y) 

n 
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(4.1) 



Where {u:,,Uy) = J^[U (x,y)] and = T[h^^ (x,y)]. But (^^, 0 ;,/) 

can be analytically evaluated [12], 



= ^[hz„] = exp 

yielding. 

0(x,y) = 

The algorithm to evaluate (4.3) is shown .schematically in 4-1, and computer 
code to perform the algorithm is given in Sections A. 1.1 and A. 1.2. 
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xP{x,y)=lFFT 

^ Final Result. , 



Figure 4-1: Hologram syntliesis algorithm 
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4.3 HoloSynth Simulator 




Figure 4-2: HoloSynth user interface 



4.3.1 Introduction 

HoloSynth assists in synthesizing liolograrns similar to those obtained from Holo- 
graphic Particle Image Velocimetry. It provides the user with a simple way to 
place any nimiber of jrarticles into a cubic volume in preparation for generating 
a synthetic hologram. The user can manually place particles anywhere in the 3- 
diniensional region as well as fill arbitrary rectangular solid regions with randomly 
distril)uted particles. 

HoloSynth uses the method for synthesizing holograms described in Section 
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□ IHoloSynth 1.0 I x 
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Figure 4-3: HoloSynth main controls 

4.2. Ill addition it prepares an input file for an alternate hologram synthesizer 
so that different synthesis algorithms may be easily compared, without modifying 
the HoloSynth source code. The output file is also useful in analyzing the results 
of reconstruction algorithms. 

HoloSynth was developed under Red Hat(R)Linux version 5.1 and tested under 
Microsoft Windows NT@. 

4.3.2 The Interface 

4.3.2. 1 The Main Controls 

Figure 4-3 shows the main controls for HoloSynth. 

Eacli button will open up a new window that controls a particular aspect of 
HoloSynth. The windows are: 

1. 3-D View of the current particle configuration. 

2. Plane View, which allows manual placement of particles on an arbitrary 
plane. 

3. Select Region, which allows selection of an arbitrary rectangular solid region 
within the view space. From this control panel you can: 

• Clear all or any part of the volume. 
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Figure 4-4: 3-D View window 

• Fill all or any part of the volume with randomly placed particles. 

4. Tlie Synthesizer Settings window whicli controls the physical characteristics 
of the simulation used for the hologram synthesis. 

Through the file menu the user can open and save particle fields, and exit 
HoloSynth. 

4.3. 2. 2 The 3-D View 

The 3-D View, shown in Figure 4-4. provides a three dimensional view of the 
current, particle configuration. There are three basic functions controlled by this 
window: the magnification (zoom) of the view, the angle and translation of the 
view, and the position of the active and two clipping planes. No direct editing 
takes place in this window, it is merely a visualization aid. 

Tli(' viewing "camera position” is controlled by the rollers and sliders to the left 
and l)elow the viewing wundow. The slider above the viewing window^ controlling 
the magnification of the camera. 

The various view control planes are controlled using the slide bars to the right of 
t he 3-D view window. They appear as colored squares drawui around the periphery 
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□ HoloSynth 1.0: Active Plane View 




Figure 4-5: HoloSynth plane editor 

of the cube. Only the active plane is visible in Figure 4-4, the clipping planes are 
hidden at the ends of their travel. On a color display the active plane is bright 
green, the rear cliiping plane is blue, and the front clipping i)lane is magenta. The 
active plane determines tire 2 -coordinate that the plane editor uses when adding 
new particles. The front and rear clipping plane restrict the range of 2 coordinates 
visible in both the plane view and the 3-D view. In all cases, particles lying on 
the active plane are visible within both views. If tire user wants to view only the 
particles on the active plane, then he places both clipping planes at the same limit 
of trav'el. 

4. 3. 2. 3 The Plane View 

The plane view projects all particles visible in the 3-D view onto a single plane. 
Particles that are on the active (green) plane in the 3-D view are highlighted green 
in the plane view. Particles can be manually positioned or deleted on the active 
l>lane. 

To add a particle, simply left click within the black portion of the active plane 
view. Accurate positioning of particles can be performed in two ways; 
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Figure 4-6: 3-D region selection 

1. Turn on the coordinate grid. At the bottom of the plane view wdndow is a 
selection list allowing presentation of grids of various densitites. High grid 
densities are com^enient at high magnifications. While a grid is displayed 
the iJarticles automatically snap to the nearest grid intersection. This can 
be controlled b}^ the Display->Snap to Grid menu item, and can also be 
toggled using the ‘g’ key. 

2. Use the real-time coordinate diplay. Any time the mouse is over the edit 
window its complete 3-D coordinate is show at the bottom of the window. 

Deletion is accomplished by dragging with the right mouse button over the 
region containing the particles to be deleted. A red selection border will appear 
showing the selected region. To delete inside the selected region, use the menu 
item. Edit->Delete Inside Selection or the delete key. Using the Edit->Delete 
Outside Selection menu item will clear all particles visible in the plane editor 
but outside the selection box. The region deleted is bounded in the x and y di- 
rections by the selection box, and in the 2 dimension by the clipping planes. If 
a i^article is visible in the plane editor, regardless of whether it is on the active 
plane, it can be deleted. 
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4. 3. 2. 4 Select Regions 

The 3-D Selection window, shown in Figure 4-6, allows the user to define a solid 
rectangular sub-volume within the 3-D view and perforin the following oirerations: 

1. Clear the entire volume. 

2. Clear all particles inside the selected sub- volume 

3. Clear all particles outside the selected sub- volume 

4. Fill the (Uitire volume with randomly placed particles. The slider at the right 
of the window controls how many particles will be added. 

5. Fill the selected sub-voluine with randomly placed particles. The slider at 
the right of the window controls how many particles will be added. 

6. Fill outside the selected sub-volume witli randomly placed particles. The 
slider at the right of the window controls how many particles will be added. 

The sub-volume is displayed as a translucent yellow box with the 3-D view. Each 
wall of the selection region is controlled by its own slider. The region will only be 
displayed in the 3-D view if the Show Selection Box button has been selected. 

4. 3. 2. 5 Hologram Synthesis 

Figure 4-7 shows the synthesizer controls. This allows easy control of the physical 
size of the simulation as well as recording plane distance and number of planes used 
to subdivide the simulation volume. The generated hologram will be displayed in a 
new window when the computations have completed. If the product is satisfactory, 
the user may save the hologram to a file using the Truevision Targa format or the 
Adol)e Photshop Raw format. 

4.3.3 Requirements 

HoloSynth is written completely in C-l— t-. The Graphical User Interface uses 
the freely available platform indei)endent user interface toolkit FLTK[13]. The 



4.3. HOLOSYNTH SIMULATOR 



45 




Figure 4-7: Synthesizer settings 



hologram synthesizer use tlie FFTW fast fourier transform library available freely 
from MIT[14], The program runs under X- Windows, and Microsoft Windows with 
OjtenGL[15] or the Mesa 3-D graphics library (an OpenGL functional equivalent) [16]. 

4.3.4 Program Architecture 

At its lieart HoloSynth is a fancy list editor, editing a list of particle coordinates 
and radii. The core data structure is. appropriately, a linked list of particle objects. 
This allows an arbitrary number of particles to be added, restricted only by the 
memory available. 

Figure 4-8 shows the schematic relationship between all program modules that 
make up HoloSynth. The complete source code to HoloSynth is presented in 
Appendix A.l. 
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Figure 4-8: HoloSynth architecture 
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Chapter 5 

Computational Reconstruction 

5.1 Onural’s Method 

In 1985, Levent Oniiral proposed a technique based on a truncated inverse filter[12, 
17], The method is summarized here, source code is presented to perform the 
algorithm in Sections A. 2.1 and A. 2. 2. 

Starting from (2.17), 



y) = U {x, y) ® h^{x, y) 

We letr U (x, y) = 1 — a (x, y) where a {x, y) is tlie opacity of the given plane. Now 
write the recording plane intensity. 



Iz {x, y) 



- i>z{x,y)il)l{x,y) 

= \{1 - a{x, y)) ® h^{x,y)f 

= 1 - a (x, y) ® hi (x, y) - a (x, y) ® (x, tj) 

+ \a{x,y)®h^ {x,y)f 
= 1 - a (x, y) ® [hi (x, y) + h^ (x, y)] 

= 1 - a(x, y) ® 2Re{/),j. (x,y)} 



(5.1) 



(5.2) 



Wh('re |a (x, y) ® hj (x, y)|" was dropped as insignificant in (5.1). Equation (5.2) 



I 
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shows the intensity modeled as a DC shifted linear system with impulse response. 



g{x, y) = 2Re{h, {x,y)} 



(5.3) 



No inverse of g (.x, y) in the normal sense exists, However, Onural show(*d in 
[12] that a series approximation to the Fourier transform of g can be found. If 
g = ^[9 (•!■, .'/)] then. 



QaI {x>s.ujy) — cos — [ul + u;^) + 






in 



k=l 



M 1 



M 



where [logo is the integer obtained by truncating the value of log 2 k and, 



Su {ujj,u)y) — ( — I)*" COS — {2k +1) {ul + tjy) 



(5.5) 



yielding. 



9M{^^y)=^ ‘ 



(5.6) 



Looking back at the hologram recording given by (5.2), 

o(x,y) ®2Re{/i, (x,y)} = \-I{x,y) 

T[o\x,y)\g {u^.LOy) = T[\-I{x,y)] 

fu{x,y) = T-^[T[l-I{x,y)]Q]^,\{oj,.u;y)] (5.7) 



yielding the original opacity function. 



5.2 Simulation Results 

5.2.1 Synthesis and Computational Reconstruction 

Figure 5-2 is a portion of a hologram generated by HoloSynth. The particle config- 
uration consisted of 2328 particles. 2000 of the particles were randomly distributed 




1 



5.2. SIMULATION RESULTS 



49 




Figure 5-1: Schematic of tlie particle configuration used to generate Figure 5-2 

tlirougliout a 5cin cube. Along the center plane parallel to the wavefront 328 par- 
ticles were positioned to spell the letters “MIT”. Figure 5-1 illustrates a schematic 
of that configuration. The synthetic hologram itself was 2048x2048 pixels, and 
the volume was subdivided into 1024 planes, for a 2 -axis discretization of 0.04mm. 
The simulated particle diameter was 25pm. The portion of the hologram shown in 
5-2 is 512x512 i>ixels and its upper left corner is 205 pixels from the left extremity 
of the original and 741 pixels from the top. Although it contains the ‘M’ from 
“MIT” , it is impossible to distinguish from the diffraction patterns produced by 
the other particles. 

Figure 5-3 through 5-8 show the results of reconstructing Figure 5-2 at various 
distances between 30.0cm and 30.5cm from the recording plane using Onural's 
method with M — 5. The images generated are the opacity functions, with black 
being transparent, white being opacpie. 

In Figure 5-3(a) the ‘M’ is clearly visible on a computer monitor with some 
sliglit clouds from out-of-focus particles. Figure 5-4 is the same hologram recon- 
structed at 30.1cm. Here there is little visible difference from the previous figure. 
This is an artifact of the way the decoder transform real valued intensities into 
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gray scale bitmap images. 

To translate from an arbitrary field of real number to a field of 8 bit integers, 
the decoder first finds the maximum and minimum values present in the real array, 
it then scales and translates each of these values to fill the range from 0 to 256. 
in order to maximize use of dynamic range available in an 8 bit gray scale image. 

The 30.1cm reconstruction has a much lower maximum intensity than the 
30.0 cm reconstruction. Volume reconstruction, where many planes will be recon- 
structed at a time, in order to reconstruct flow fields, will have to keep track of 
this scaling issue between planes. 

The important thing to note about this series of decodes, is that while the 
source hologram was constructed using a 2 discretization of 0.04mm, the recon- 
struction algorithm 2 -discrimination was of the order of millimeters. Thus, the 
original image could have been discretized into 64 steps, or 0.78mm slabs, requiring 
l/16th the computation time with little or no loss of reconstruction accuracy. 

5.2.2 Image Postprocessing 

The next step in reconstruction is to automatically detect particles in each plane 
and register their x, y, and 2 coordinates, rather than simply displaying an im- 
age. Repeatedly reconstructing at different depths and detecting the particles will 
replace the time consuming 3-D scanning. 

In Figures 5-3 through 5-8 subfigure (b) and (c) show the results of applying 
.some simple image processing techniques to the original reconstructions to enhance 
the output. Subhgure (b) in each figure applies edge detection', which simply 
Hlters the image marking any region where the gradient of the image exceeds 
a certain threshhold. In the case of the grayscale images being used here, the 
gradient, is simply the difference in grayscale value between adjacent pixels. [18] 
Subfigure (c) in each Hgure ‘threshholds’ the image, changing a pixel to white 
(value 255) or black (value 0) depending on whether or not the i)ixel value exceeds 
a preset threshhold.* 



'The image processing was accomplished with the freely available GNU Image Manipulation 



5.2. SIMULATION RESULTS 



51 



The combination of edge detection and threshholding can be used very effec- 
tively to automatically detect particle i)ositions. At the end of the image process- 
ing, there are only regions of black (value 0) with very small regions of white (value 
255). The white regions correspond to particles and can easily be automatically 
tabulated and mapped into x and y coordinates. 



Progiaiii (GIMP, http://www.gimp.org), which is very similar to Adobe Photoshop(fi). 



5.2. 
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Figure 5-2: Portion of a synthetic hologram. The original synthetic hologram wius 
2048x2048 pixels. This portion is 512x512 pixels, cropped from the orig- 
inal hologram to make the images j)resentable in printed form. 
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(a) Raw Reconstructed Image (b) Image after enhancing with edge de- 

tection. 




(c) Image after enhancing with edge de- 
tection and threshholding. 



Figure 5-3: Figure 5-2 decoded at a depth of 0.300in. 
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(a) Raw Reconstructed Image (b) Image after enhancing with edge de- 

tection. 




(c) Image after enhancing with edge de- 
tection and threshholding. 



Figure 5-4: Figure 5-2 decoded at a depth of 0.301ni. 
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(a) Raw Reconstructed Image (b) Image after enhancing witli edge de- 

tection. 




(c) Image after enhancing with edge de- 
tection and threshholding. 



Figure 5-5: Figure 5-2 decoded at a depth of 0.302m. 
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(a) Raw Reconstructed Image (b) Image after enhancing with edge de- 

tection. 




(c) Image after enhancing with edge de- 
tection and threshholdmg. 



Figure 5-6; Figure 5-2 decoded at a depth of 0.303m. 
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(a) Raw Reconstructed Image (b) Image after enhancing with edge de- 

tection. 




(c) Image after enhancing witli edge de- 
tection and threshholding. 



Figure 5-7: Figure 5-2 decoded at a depth of 0.304in. 
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(a) Raw Reconstructed Image 



(b) Image after enhancing with edge de- 
tection. 




(c) Image after enhancing with edge de- 
tection and threshholding. 



Figure 5-8: Figure 5-2 decoded at a depth of 0.305ni. 
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Chapter 6 

Conclusions and 
Recommendations 



It is clear that computational reconstruction of HPIV is practicable. The proce- 
dure that Onural proposed in 1985 was not then practical due to the computational 
jrower commonly available. Modern computing can easily handle realistic holo- 
grams and tlie numeric techniques of normal PIV could be extended to deal with 
the particle identification problem in HPIV. The simulator presented in this thesis 
will be useful in any setting when attempting to validate particular reconstruction 
algorithms. 

Future work should concentrate on extending the reconstruction prosses fully 
into three dimensions and developing particle identification techniques. Automat- 
ically determining the velocity vector for a particle is the ultimate goal. 
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Appendix A 



Source code 



A.l HoloSynth: Hologram Synthesizer. 

This is the source code for HoloSynth, a program for easily creating particle configuration 
and holograms for use calibrating computational reconstruction algorithms. 

A. 1.1 OnuralSynthesizer.hpp 

#ifndef ONURALSYNTHESIZER_H 

/ 1! include guard 

#def me ONURALSYNTHESIZER_H 1 

#include ’’Synthesizer .hpp” 

#include ’’ParticleList .hpp” 

#ifndef FLMATH 
#define FLMATH 
#include <FL/math.h> 

#endif 

#include RGB Image .hpp’* 

/^!\ class OnuralSynthesizer 
* 

^\brtcf The algorithm use for actually generating the hologram 
* 

* The forward diffraction algorithm implement a method described in 

* ''Digital Decoding of In-Line Holograms^^ a PhD thesis by Leven Onural 

class OnuralSvnthesizeripublic Synthesizer 

{ 

public: 

/// created with a reference to the current particle list 
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OiiuralSyntlu\sizer(ParticleList *pl); 

~OiiuralSynthesizer(); 

j j\ Generate, organized all the FFTs and diffraction calculations. 
int generate(); 

private: 

RGBIinage * hologram ; 

void diffract (fftw_coinplex *image); 
void clearField(fftw_coinplex *ffeld, fftw^real val); 
void plotParticle(particle *c.fftw_complex *plane); 
void plot(int p_X, int p.Y, fftw.complex *plaiie); 
void addto(ff'tw_complex *al, fftw.complex *a2); 
void MagtoRGB( fftw.complex *data); 

void cleanupO; 

fftwnd.plan fftwFwdPlan; 
fftwiid_plan fftwBwdPlan; 

fftw.complex *oplane; 
fftw.complex *hplane; 

}; 



#endif 



A. 1.2 OnuralSynthesizer.hpp 

#include <FL/F1.H> 

#include ’’ProgressUI .hpp” 

#include "ScrollImageUI . hpp" 

#include <FL/Fl_Image . H> 

#include ” OnuralSynthes izer . hpp” 

#include ’’StatusBox .hpp’’ 

#def ine MIN(a, b) (((a) < (b)) ? (a) : (b)) 
#def ine MAX(a, b) (((a) > (b)) ? (a) : (b)) 



OnuralSynthesizer::OniiralSynthesizer(ParticleList *pl) 
.•Synthesizer (pi) 

{ 

} 

OmiralSynthesizer::^OnuralSynthesizer() 

{ 



cleamipO; 




i 
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} 

int OimralSviitliesizer::geiierate() 

{ 

int progresscount=0; 

progressUI— >setMax(5-f 5*d_Plcines); 
progressUI — >upclate( 0) ; 
progressUI — > reset ( ) ; 

{)rogressUI -^show( ) ; 

Fl::wait(0.0); 

particle *curr; 

int pixels, il, pcount, Ncalc, Mcalc; 
fftw.real zgrain, Zinin, Zinax; 



* Initialization. 

3^cs|c5|c5^c5tc5|C5^cic3|e5fc:4::+:5^c5f::+:5ic5ic:+:5tc5^(:5|<:5f:5|e5^c5^c5^c5|e5f:^5^<::+:5ic:+C3fc:fC5|c5^c3f:5f::+::+:3^c:4c:l|c:+::4c3tc*:+^:+:*:+^ 



pixels=2*(2*N/2+l)*2*M; 

Ncalc=2*N; f j Double the size for padding out the edge effects. 
Mcalc=2*M; 



if(progressUI— >setStatus("Allocating memory. . .")){ 
cleanupO; 
progressUI — > hide ( ) ; 

return 0; 

}: 

hplane=(fftw_complex *) calloc( pixels, sizeof(fftw_complex)); 
oplane=(fftw_complex *) malloc( pixels*sizeof(fftw_complex)); 

progressUI— >update(++progresscount); 

zgrain=(fftw_real)l/(2*d_Planes);//T/ie thickness of a HALF plane 

if(progressUI— >setStatus("Generating FFTW Plans ... ")){ 
cleanupO; 
progressUI— >liide(); 
return 0; 

}; 



fftwFw(lPlan=fftw2d.create.plan(Ncalc, Mcalc, FFTW_FORWARD. 



A.l. HOLOSYNTH: HOLOGRAM SYNTHESIZER 



63 



FFTWJNJ^LACE); 

fftwBwdPlan=fitw2cLcreate-plaii(Ncalc, Mcak\ FFTW_BACKWARD, 
FFTWJNJ^LACE); 

progr essUI —> update (+ + progr esscoimt ) ; 

1 1 Rememher we will eventually need to split out the vector end of things. 
pList^sortListOnZ(); 

curr=pList^first(): 

for(il = l;il<=d_Planes;il++){ 



Fl::wait(0); 

if(progressUI— >setStatus("Clearing Object Plane... ")){ 
cleanupO; 
progressUI ^ hide( ) ; 

return 0; 

}; 



clearField(oplane, 1.0); 

progressUI— >update(++progresscount); 

Zmin=(iTl)*2*zgrain-0.5; 

Z=(Zmin+zgrain)*depth;/ /Physical Z 
Zniax=Zmin-)-2*zgrain; 

pcount=0; 

if(progressUI— >setStatus("Plotting Particles. . . ")){ 
cleanupO; 
progressUI ^ hide ( ) ; 

return 0; 

while((curr!=NULL)&:&:(curr— >Z>Zmin)&:&:(curr— >Z<=Zmax)){ 
plotParticle(curr,oplane); 
pcount++; 
curr=curr^next; 

} 

progressUI ^ update( + +progresscount ) : 
if(pcount>0){ 



if(progressUI->setStatus("Perf orming Forward FFT...")){ 
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}; 



rleaiiupO; 

[) r ogress U I ^ li i d e ( ) ; 

return 0; 



fftwnd_oiie(fftwFwdPlan, oplane, 0); 

I^r ogr essUI ^ u p d at e ( + -f pr ogresscount ) ; 

if(progressUI^setStatus(”Dif f racting. . . ")){ 
(deanupO; 
progressUI— >hide( ) ; 

return (); 

}; 

diffract (oplane); 

progressUI— >update(++progresscount); 

if(progressUl^setStatus( "Adding. . .")){ 
cleanupO; 
progressUI^hide(); 

return 0; 

}; 

addto(hplane,oplane); 

progressUI— ^update(+-f progresscount); 

} 

progresscount=5*il+2; 

progressUI^update(p)rogresscount); 

} 

hplane[0].re=10000()0; 

progressUI—)>setStatus(" Performing Backward FFT..."); 
progressUI— >update(++progresscount); 
fftwnd_one(fftwBwdPlan, hplane, 0); 
progressUI—>setStatus( "Calculating Intensity. . • 

MagtoRGB (hplane); 

f I The hologram is finished, now display it... 

FlJmage *image=new FUmage(holograni^data(), N, M); 
progressUI—>setStatus( "Displaying. . ."); 
progressUI^npdate(++progresscount); 

ScrolllmageUI *si=new ScrollImageUR image , 0,0, 500,500); 
si— )>sliow(); 

progressUI ^ li ide( ) ; 
cleanupO; 



} 



return 1; 



i 
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void OmiralSynthesizer;:addto(fitw -Complex *a, fftw_coini)lex *b) 
int il; 

for (il=0: il < ( 4 * N * M ); il++){ 

if((fabs(b[il].re)>=0.0000001)&&(fabs(b[il].im)>=0.0000001)) 

{ 

a[il].re -h= b[il].re; 
a[il].im -f= b[il].im; 

} 

} 

} 

void OiiuralSynthesizer::clearField(fftw_complex *field, fftw_real val) 

{ 



int il; 

for(il=0; il < ( 4 * N * M ); il-{-H-){ 
field [il]. re = val; 
field[il].im = 0.0; 



} 

void OnuralSyiitliesizer::plotParticle(particle *c,fftw_complex *plane) 

{ 

int X, y, r; 



x=(int)((c^X + 0.5)*N); 
y-(int)((0.5 - c-^Y)*M); 
r=(int)(c->R*N); 



sv^itch (r) 

{ 

case 6: 

plot(x-l, y+1, plane): 
plot(x-l, y-1, plane); 
plot(x+l, y-1, plane); 
plot(x+l, y+1, plane); 
case 2: 
case 3: 
case 4: 
case 5: 



plot(x, y-1, plane); 
plot(x, y-hl, plane); 
plot(x-l, y, plane); 
plot(x-{-l, y, plane); 
default: 

plot(x, y, plane); 
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} 



} 

void OmiralSyntliesizer;:plot(int p_X, int p_Y, fftw.complex *plaiie) 

{ 

int pos=2 * p-Y * N + p_X; 
if(pos>0){ 

p 1 an e [ p os] . r e = 0 . 0 ; 

])lane[ pos].iiii=().0; 

} 



} 

#define F.PI 4.0*3.1415926535 

void OnuralSynthesizer::diffract(fftw .complex *irnage) 

{ 

int X, y, Ncalc, Mcalc; 
long ik; 
fftw.real Iz, cc; 
fftw.real wx, wy, a, ca, sa; 

Ncalc=2*N; 

Mcalc=2*M; 

lz=L*(distance-Z); 

ik=-0: 



/* FFTW puts frequency space origin at the upper right corner of the 

* army. Positive frequencies build from the upper left, corner 

* diagonally down through the upper left quadrant. The other corners 

* of the array the array are ALSO at the origin, so the very center 

* of the array is the HIGHEST frequency. 

* 

* In order to use the diffraction filters we MUST take that into 

* account, that is what the strange trigraph statements take care of 

* below */ 

for (y=0;y<Mcalc;y4-+) 

{ 

wy = (y<M) ? y/height : (y-Mcalc)/height; 

// wy — (y<M) ? y/height : (Mcalc-y) /height; 
wy * = wy; 

for(x=0;x<Ncalc;x+-f ){ 

wx=(x<N) ? x/width : (x-Ncalc)/width; 

II w:r^(x<N) ? x/width : (Ncalc-x) /width; 

a=lz*(wy'fwx*wx) *3. 141 5926535; 

/ / a.= Iz^ ( wy-\- w:r^ wx)/F, PI; 
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ca=cos(a); 

sa=sin(a); 

cc=irnage[ik].re*ca-iniage[ik].iiii*sa: 

image[ik].iiii=image[ik].re*sa+image[ik].iiii*ca; 

image[ik].re=cc; 

ik++; 

} 

} 

} 

void OriiiralSyntliesizer::MagtoRGB(fftw_complex *data) 

{ 

int i. j, ij, grey, Neale, Mcale: 
fftw.real v, viiiax, vmiii, vl. v2: 



Ncalc=2*N; 

Meale=2*M; 



hologram=new RGBImage(N,M); 

vmax=0.(); 

vmiii=lcG; 

ij=0; 

/ j Scale the result to take maximum advantage of 8bit range 
for(i=0:i<Meale*Ncalc;i+H-){ 
vl=data[i].re; 
v2=data[i].im; 
v=sqrt(vl*vl+v2*v2); 

vmax=MAX(v, vmax); 
vmin=MIN(v, viiiin); 

} 

vmax=l/ (vniax“vmin); 



RGBColor c; 
for(i=0;i<M;i++){ 
for(j=0;j<N;j + + ){ 
ij=i*Nealc+j; 



vl=data[ij].im; 

v2=data[ij].re; 

v=(sqrt(vl*vl+v2*v2)-vmin)*vmax; 

data[ij].rc=v; 
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} 



} 



grey=(int)ttoor((f{t.w.real)(v*255.0)): 
c.B=grey; 
c.R— grey; 
e.G=grey; 

liologram— ^plot(j, i, c); 



} 

void OiiuralSynthesizer::cleanup() 

{ 

free(hplaiie); 

free(oplane); 

fftwnd_destroy_plan(fft\vFwdPlan); 

fftwiKLdestroy>plaii(fftwBwdPlan); 

} 




1 1 > 



A.2. DECODER; ONURAL RECONSTRUCTION ALGORITHM 



69 



A. 2 Decoder; Onural Reconstruction Algorithm. 



This is the source code for an Onural reconstruction algorithm. 



A. 2.1 OnuralDecoder.hpp 



#ifndef ONURALDECODER.HPP 
#define ONURALDECODER.HPP 1 



#include ’’Decoder .hpp'" 

#include ’’ProgressUI .hpp" 

#include <fftw.h> 

class OnuralDecoder: public Decoder 



{ 



public: 

OnuralDecoder 0 ; 

virtual uchar* decode(uchar *data, int N, int M, 

float W, float H, float L, 
float D); 



private: 

int dJmageN; 
int d_CalcN; 
int dJmageM; 
int d-CalcM; 

float d JmageW; 
float dJmageH; 
float d JmageLambda: 
float d_decodeDist; 

fftw_complex *d_Field; 

uchar * d_data; 



void OnuralDecoder::noDClquad(fftw_complex *field, fftw_real average Value); 

fftw.real OnuralDecoder: :aprinv(fftw_real zz, int Mstage); 

void OnuralDecoder::onural(fftw .complex *hg, 

int N, int M, 

fftw_real L, fftw.real Z, 

fftw_real X, fftw.real Ydnt Mstage); 



fftw.real Magnitude(fftw_complex *field): 



void uchar2field(); 
uchar *tield2uchar(): 
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fitAv.real *cli{)Field(); 



iichar *OnuralDecoder::niaiiiLoop(); 
ProgressUI ^progress; 

}: 

#endif 



A. 2. 2 OnuralDecoder.cpp 

#include "OnuralDecoder .hpp” 

#include <fftw.h> 

#include ’’StatusBox .hpp” 

#include <FL/F1.H> 

#include <math.h> 

#def ine PI 3.141592G5359 

#def ine TPl 6.28318530718 

#define MIN(a, b) (((a) < (b)) ? (a) : (b)) 

#def ine MAX(a, b) (((a) > (b)) ? (a) : (b)) 

OnuralDecoder::OnuralDecoder() 

{ 

d JmageN=0; 
d JinageM=0; 

d JmageW=0: 
d JinageH=0; 
d J in ageL amb d a= 0 ; 
cLdecodeDist— 0; 
dMata=0; 

progress=new ProgressUI(); 

} 

ucliar *OnuralDecoder::decode(uchar* data, int N, int M, float W, float float L. 
float D) 

{ 

d JiiiageN=N; 
d_CalcN=2*N; 

dJinageM=M: 

d_CalcM=2*M; 

dJniageW=\V; 

dJmageH=H; 

dJmageLambda=L; 

d_decodeDist=D; 

d_data=data: 
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return maiiiLoop(); 



/* Remove the DC component from only the first Quadrant, the other 
* three quadra, nts are ALREADY ZERO^ j 

void OiuiralDecoder::iioDClquad(fftw_complex *field, fFtw_real average Value) 
{ , 

int i, j, pos: 

for(i=0; i<d_ImageM: i++){ 

for(j=0; j<dJmageN; j++){ 
pos=i*d_CalcN +j : 
field [pos]. re -= average Value; 

} 

} 



fftw_real OniiralDecoder::aprinv(fftw_real zz, int Mstage) 

{ 

fftw.real cc,ss,powKa; 
int k,powK: 

cc=cos(zz); 

ss=(); 

if(Mstage>l){ 
ss=2*sin(zz); 
cc+=ss*sin(2*zz); 
ss* =2*cos(2*zz); 

} 

for (k=3;k<=Mstage;k++){ 

powK=4^(k-3); // powK=2^*(k-l)^ fine for small k, and FAST 

powKa=(fft\v_real)powK*zz; 

cc+=ss*sin(powKa); 
ss* =2*cos(powKa); 

} 

return 2*cc/Mstage; 



} 

void OiiiiralDecoder::oiiural(fftw_complex *lig, 
int N, int M, 
fftw.real L, fftw^real Z, 
fft\v_real X, fftw.real Y,int Mstage) 



{ 
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int X, y. Ill 12. iil2; 
long ik; 

fftw.real alpha. Iz. wx2, wy2; 

lz=L*Z: 

ml2=M/2; 

nl2=N/2: 



ik=(); 



for (v=0;v<M:v++)//M 

wy2 = (y<ml2) ? y/Y : {y-M)/Y; 
wy2 * = wy2: 
for(x=0;x<N;x++)//A^ 

{ 

wx2 = (x<nl2) ? x/X : (x-N)/X; 
wx2 * = wx2; 

alpha=PI*lz*(wy2+wx2); 

lig[ik].re* =aprinv(alpha, Mstage); 
ikH — h; 

} 

if( progress— >update(3) ) { 
progress— )*hide(); 

return; 



} 

} 



fftw.real OniiralDecoder::Magnitude(fftw.complex *field) 

{ 

int i: 

fftw.real v,suinV, vl, v2; 

siiinV=0: 

for(i=(): i < d_CalcN*d.CalcM; i++){ 
vl=field[i].re; 
v2=ffeld[i].im; 
v=sqrt(vl*vl+v2*v2); 
lield[i].re=v; 
field[i].im=0; 
siimV+=v: 
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} 



return sumV/(d_CalcN * d.CalcM); 



} 

void OmiralDecoder::iichar2tield() 

{ 

int arraysize=d_CalcN*d_CalcM*sizeof(fftw_coinplex); 
int fieldpos, iicliarpos, x, y; 
iidiar niinval, maxval; 
fitw_real scale; 

d_Field=(fftw_coin})lex *)calloc(arraysize. 1);/ / calloc clears to zero... 

iiiinval=0; 

maxval =255; 



for(y=0; y< dJmageM; y-f + ){ 
for(x=0; x<dJmageN; x++){ 

ucharpos=3*(x + y*dJmageM); 
minval=MIN(minval, d_data[ucharpos]); 
maxval=MAX(maxval, d-data[ucharpos]); 

} 

} 

scale= 1.0/ (fftw.real) (inaxval-minval) ; 

for( y=0; y< dJmageM; y++ ){ 
for( x=0; x<dJmageN; x-f+){ 

fieldpos=x + y*d_CalcN; 
ucharpos=3*(x -f y*dJmageN); 

d_Field [fieldpos]. re= ((fftw.real) (d_data[ucharpos]-minval))*scale; 

} 



} 

iicliar *OnuralDecoder::field2iichar() 

{ 

int arraysize=3*dJmageN*d_ImageM*sizeof(uchar): 
int fieldpos, ucliarpos, x, y; 
fftw.real miiival, maxval, \% vre, vim; 
fftw.real scale, val; 

iichar *retval; 

retval=(iichar* )calloc(arraysize, l);//calloc clears to zero... 
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iniiival=le20: 

inaxval=0; 

for(y=(); y< clJmageM; y++ ){ 
for(x=0; x<dJmageN; x++){ 
fieldpos=x + y*d_CalcM; 
vre=d -Field [fieldposj.re: 
vim=dJField[fieldpos].im; 
\~sqrt(vre*vre-j-vim*viin); 

1 n i n val = M I N ( n i i rival , v ) ; 
maxval=MAX(maxval, v); 
d_Field [fieldpos].re=v; 



} 

} 

scale=255.0/ (inaxval-minval); 

for(y=0; y< dJmageM; y++ ){ 
for(x=0; x<dJmageN; x++){ 
fieldpos=x + y*d_CalcN; 
iicharpos=3*(x + y*dJmageN); 

val=(uchar)((d_Field[fieldpos].re-minval)*scale); 
retval[ucharpos] = (uchar)val; 
retval[ucharpos+l] = (uchar)val; 
retval[ucharpos+2] = (uchar)val; 



} 

} 

return retval; 



iK'har* OiiuralDecoder::mainLoop() 



tftw_real DC; 

int progresscount—0; 

I ) r og r ess -4 set M ax ( 5 ) ; 
j)rogress-4update(0); 
j)i ogress^ reset 0; 
j)rogress-4show(); 

if(progress—)>setStatiis( "Transforming Input. . .")){ 
progress-4liide(); 
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return 0; 



fitwnd_plan ForwardFFTPlan, BackwardFFTPlaii; 
iicliar2field(); 

i:>r ogress ^ ii p date ( -h -h j:>r ogresscoi iiit. ) ; 
if(progress^setStatus( “Making FFTW Plans... “)){ 
j)rogress— >hide() ; 

return 0: 



ForwardFFTPlan=fftw2d_create_plan(d_CalcN, d_CalcM, 

FFTWJORWARD, 

FFTW_ESTIMATE|FFTWJN_PLACE); 



BackwardFFTPlan=fftw2d_create^plan(d_CalcN, d.CalcM, 

FFTWJ3ACKWARD, 

FFTW_ESTIMATE|FFTWJN_PLACE); 



DC=4*Magnitiide(d-Field); ffThe 4 handles the doubling of arraysize 



iioDClquad(d-Field, DC); 

progress-^iipdate(-h-hprogresscount); 

if(progress— >setStatus("Calculating Forward transform. 

progress->hide(); 

return 0: 

}; 



fftwnd_one(ForwardFFTPlan, dJield, 0); 

d_Field[0].re=0; 

d_Field[0].im=0; 

})rogress^iipdate(H-+progresscoiint); 
if(j)rogress->setStatus(“Starting Onural reconstruction. . .“)){ 
progress-^ liide(); 
return 0; 

}; 

oniiral(d_Field, d_CalcN, d.CalcM, d JmageLambda, d.decodeDist, 
dJmageW, dJmageH, 5); 
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progress— > update( + + progresscoimt ) ; 
if(progress—>setStatus( "Calculating the Reverse Transform..."))! 
progress—^ liide(): 

return 0: 

}; 



fitwii(Lone(BackwardFFTPlan, d_Field, 0); 

progress-^ update (++ progresscount ) ; 
if(progress->setStatus( "Complete"))! 
progress— > hide( ) ; 

return 0; 

uehar* decoded=field2uchar(); 

free(d_Field); 

progress->hide(); 
return decoded; 
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